New theoretical techniques are being developed and characterized. These efforts are usually coupled with software development, and involve the systematic testing and evaluation of new ideas. Ongoing development of the Action-CSA method We developed a new computational approach to find multiple conformational transition pathways of biological macromolecules with fixed initial and final states via global optimization of Onsager-Machlup (OM) action. For efficient global optimization of OM action, we used the conformational space annealing (CSA) algorithm and the modified classical action. This approach successfully samples not only the most dominant pathway but also other possible ones without initial guesses of reaction pathways. It was demonstrated that the rank order of actions and transition time distribution of multiple pathways identified by the new approach are in good agreement with those of molecular dynamics simulations. The new method successfully finds multiple folding pathways of the small protein FDS-1. The lowest action folding pathway obtained with Action-CSA shows that the helix of FSD-1 folds earlier than the -sheet part. This pathway is consistent with the experimentally identified folding mechanism of FSD-1. In addition, we are currently developing an algorithm to perform direct optimization of OM action using the Hessian matrix of a potential. Currently, we are investigating pathways of alpha-beta transition caused by nine-residue small peptide whose sequence is YQNPDGSQA. Benchmark of polarizable force field models using QM/MM polarization energy We compare different polarizable force field models for estimating the polarization energy in hybrid QM/MM calculations. Two models (PAC and PAD) have been formulated based on QM response kernels, and closely resemble the fluctuating charge polarizable model and the induced dipole polarizable model respectively. An empricially parametrized force field, the Drude polarizable force field, is also tested. Our results suggests that polarization effects, including both local charge distortion and intramolecular charge transfer, can be well captured by induced dipole type models with proper parametrization. Development of the MPID polarizable force field We developed a new polarizable force field based on mapping the electrostatic model optimized in the context of the Drude force field onto a multipole and induced dipole (MPID) model. Condensed phase simulations on water and 15 small model compounds show that without any reparametrization, the MPID model yields properties similar to the Drude force field with both models yielding satisfactory reproduction of a range of experimental data. We subsequently show that the MPID model also works for biomacromolecules including proteins and nucleic acids. With the new MPID force field, more than 15 years of development of the Drude polarizable force field can now be used without the need for dual-thermostat integrators nor self-consistent iterations. Constraint free energy calculation A good way to increase the accuracy and precision of free energy calculations is to use constraints and to calculate the free energy costs of constraints during post-processing. A new functionality in CHARMM has been implemented to compute the free energy when applying or removing constraints on arbitrary degree of freedom, which will be done as additional explicit steps in the free energy cycle. With this constraints for free energy methods employed, the phase space overlap between ensembles can be highly increased, which is required for accuracy and convergence. The new techniques focus on hard degrees of freedom and use both gradients and Hessian Estimation. As Hessian matrix computation is the most time-consuming step, we introduced an approximation which saves much machine time. As a byproduct for debugging this functionality, we also implemented an optimizer for the geometry of any specific set of molecular fragments. An Eighth-Shell(ES) implementation of P21 periodic boundary condition in CHARMM Spatial domain decomposition methods, for evaluating pairwise particle interactions, function by assigning a specific region of space to each processor. For moderate sized clusters (with 150 processors), eighth-shell method has been shown to provide the best performance. The parallel engine of CHARMM, DOMDEC, however, implements only the P1 periodic boundary condition. In this work, we extend the implementation to include P21 periodic boundary condition. P21 space group is represented by (-X,Y+,-Z) i.e. it has a half-screw symmetry. We modify the import region of DOMDEC while keeping the import volume unchanged and hence obtain parallelization efficiency similar to Eighth-Shell method. This method is useful for membrane simulations as it allows the movement of lipids between the bilayers, thus balancing the chemical potential between them. Free energy calculations using neural-network based potentials The extremely high computational cost of QM energy calculations limits their use on large biological systems of interest. Semi-empirical methods and classical MM force fields trade accuracy for speed and can give approximate results. However, deep learning of QM potential using neural networks has recently been shown to have similar accuracy to DFT while taking up to six orders of magnitude lesser time. In this work, we are developing a method to perform free energy calculations using neural networks based potential. We represent the molecule using vectors describing the radial and angular environment. An additional per-atom lambda variable scales the contribution of the atom to the evaluated potential energy. Lambda values of 0 and 1 represent the end states with intermediary values representing the alchemical path of transformation. The method is being implemented in OpenMM and uses Tensorflow to learn the model. Improving sampling in constant pH simulations through use of reservoirs Determining accurate pKa values with constant pH simulations requires extensive conformational sampling. At this point constant pH simulations are still much slower than regular MD simulations with fixed charges. We have thus developed a method to use reservoirs of structures generated with fixed charge MD simulations to improve sampling in constant pH simulations. The reservoir structures are then fed into the constant pH simulations in a statistically correct way. We have tested this method on small systems and show excellent agreement in calculated pKa values and sampled conformations between the new method and the standard constant pH method. For larger systems where conformational sampling is an issue, this method should greatly increase the accuracy of pKa calculations. Other ongoing efforts include: A virtual mixture simulation approach for constant pH simulation in explicit water Improving sampling in constant pH simulations through use of reservoirs A generalized self-guided Langevin dynamics simulation method